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1. INTRODUCTION 


The report presents an introductory discussion of the mathematics 
pertaining to the directional analysis of ocean waves. The presenta- 
tion is tutorial in form but does require a reasonably complete mathe- 
matical background; a background equivalent to that required in reading 
Kinsman's textbook Wind Waves (1965). 


The level of the presentation is moderate at the beginning. The 
level picks up rapidly toward the middle but there should be sufficient 
detail and redundancy in the mathematics to allow the reader to follow 
the development without having to rediscover too many omitted steps. 

It is in this sense that the report is tutorial. In some places the 
mathematical development is intuitive rather than rigorous. This is 
deliberate in order to provide insight and understanding. In most such 
cases, references to rigorous reports are given. 


The development is reasonably detailed so that the interested 
reader may apply the methods presented and use the report as an entry 
point into the rigorous theory of the directional analysis of ocean 
waves. In this respect, if the report serves as a bridge across the 
gap between a handbook and a rigorous and sparse theory on the subject 
then the objective of the report will have been fulfilled. 


The report first presents an intuitive development of a sea sur- 
face model that assumes the sea surface to be a two-dimensional random 
process definable in terms of a directional power spectrum. A discus- 
sion of the space and time covariance function and its relationship to 
the directional power spectrum follows. Both one- and two-sided power 
spectra are discussed; however, the main development is in terms of 
the two-sided spectrum. Next, the relationship between the power and 
cross power spectrum for two fixed locations and the sea surface 
directional spectrum is developed. Explicit relationships for the 
special cases of an isotropic sea and a single wave of a given direc- 
tion and frequency are then obtained. The related topic of the direc- 
tional resolving power of an array of wave transducers is then 
presented. 


Using the preliminary developments as a basis, several methods for 
the directional analysis of ocean waves based on the information 
obtainable from an array of wave transducers are presented. The methods 


are basically a direction finder technique, a least square single-wave 
train fit, and a Fourier-Bessel expansion fit. In conclusion, a 
generalized Fourier expansion method is suggested. Extensive results 
of the application of the least square single-wave train fit are pre- 
sented in Appendix A. Appendix B is a FORTRAN II listing of a program 
for this analysis. 


2. WAVE MODELS 


In its simplest form an ocean wave can be thought of as a single 
frequency, sinusoidal, infinitely long crested wave of length \, moving 
in time over the ocean surface from a given direction 6. Such a wave 
is illustrated in Figure l. 


Spatial Frequency N(usv, ty) Spatial Frequency Along 
Along v axis is zero. y axis is m = K sin 8 


Direction of @ Wave 
Travel in Time t 


u Spatial Frequency 
along u axis K = 1/r 


Spatial Frequency 
along x axis is 
£=K cos @ 


FIGURE 1. SIMPLE OCEAN WAVE 


Assume that the wave is frozen in time over the surface (the Xyy 
spacial plane). The coordinates (u,v) are a 9 degree rotation of the 
(x,y) coordinates. The positive u axis lies along the direction from 
which the wave is traveling. The wave surface n(u,v), shown frozen in 
time in Figure 1 can be described mathematically by 


n(u,v) = cos(2nKu + 279) (2a) 


where K = 1/\ is the wave number of spacial frequency in cycles per 
unit length along the u axis, and 27d is a spacial phase shift. 


To make the wave move in time across the spacial plane with a 
time frequency f = 1/p, where p is the wave period, it is necessary to 
add a time part to the argument of the cosine function in the model 
above. The time part is a phase shift dependent only upon time. As 
time passes, the time part changes causing the cosine wave to move 
across the (u,v) plane, in this case the ocean surface. Adding the time 
part we get (where 27) is a fixed time phase shift) 


n(u,v,t) = cos(2nKu + 276 + 2n ft + 2my). 
If we combine the effect of the ¢ and wt phase shifts as a = $ + yy, 
we get 

n(u,v,t) = cos(2n(Ku + ft + a)) (252) 
as a simple model of a sinusoidal wave moving in time over the ocean 


surface. 


Since the coordinates (u,v) are a rotation of the coordinates (x,y) 
through an angle of 6 degrees, we know 

ul =x cos 8) y san) ¢ 

V = —-x sin 0 + y ‘cos 0: 
Using the above relations, and letting 2 = k cos 6 and m= K sin 6 


be the spacial frequencies along the x and y axes, respectively, we 
have 


n(x,y,t) = A cos(2n(2£x + my + ft + a)) @23)) 
as a model for a wave of height 2A moving from a direction 


6 = arctan (m/2) 


with a phase shift of 2ma. A wave crest of such a wave system is infi- 
nite in length. A crest occurs at a set of points (x,y,t) which satisfy 
the relation 


gx + my + ft = a constant = (n - a) 


where n = QO, -1, tl, -2, +2, ... . Each value of the index n relates to 
a particular crest. The intersections of the crests with the x and y 
axes move along the respective axes with time velocities V, = -f/2 

and Vy = -f/m. This follows from the differential expressions 


D,@)= D, [aes _ mo -it] -t (2.4) 


x M-% = LH t 
DweD[ Ms te ft J.-£ es 


obtained from the wave crest relationship given above. 


From Euler's equation we know that cos y = (Exp(iy) + Exp(-iy))/2. 
If we consider y as 2n(2£x + my + ft + a) we can write 


n(x,y,t) = 1/2A Exp(i2n(2x + my + ft + a) + 
1/2A Exp (i2n(-2£x - my - ft - a)) @z6) 


In the above we have introduced the notion of negative time frequencies. 
This makes it possible to express an elementary wave in the mathemati- 
cally convenient form 


n(x,y,t) = a Exp (i27(2x + my + ft + a)) (2.7) 


where a = 1/2A. In the real world a complex wave of this type implies 
the existence of another wave n¥*(x,y,t) which is the complex conjugate 
of n(x,y,t) above. This complex conjugate is given by 


n*(x,y,t) = a Exp (-i2n(2£x + my + ft + a)) 
=a Exp (i2n[(-2)x + (-m)y + (-f)t + (-a)]). (2.8) 


The fact that negative frequencies are considered is explicit in the 
above relation. 


A property of the above model, which will be used later in connec- 
tion with the directional analysis of waves from measurements obtained 
from an array of detectors, is expressed by the equation for the phase 
difference of two measurements made at two different points in space 
and time. Assume we know the value of n(x,y,t) at the three-dimensional 
coordinates (Xo,yo,t,.) and (xo+X, yotY, t +I), where X, Y, and T 
are constants. The phases at the two points are given by 


(x5 »¥52t,) = LX, ate Wb ate ft, + a, (2.9) 
(x, + X, Vg + Y, ty + T) = (x, + X) + m(y, + Y) + F(t) +T) +a 


(2.10) 
This gives a phase difference of 
Ad = (2X + mY + fT). @ way) 


To obtain a more complicated wave system consisting of many waves 
of various frequencies and directions, we can linearly superimpose (add 
up) many waves of the form given above. If we do this, we can write 


N 
1 (x.9.t) =>, Exp (i 20 (InX+ nd +f, ¢ +05), (2ea2) 
NB | 


For this wave system to be real, the terms must occur in complex con- 
jugate pairs as indicated above. 


For completeness, consider a model for an infinite but countable 
number of distinct (discrete) waves and write 


n(x,y,t) — S a, Exp (i an (I, K+ mY +f,t + %y)) : (DNs) 
ne \ 


Again the terms must occur in complex conjugate pairs for the wave 
system to be real. This will be assumed to be the case in future 
discussions. 


A model for a wave system in the case where energy exists for con- 
tinuous intervals of frequency and direction should be considered. In 
particular, consider the general case of continuous direction from 
0 to 27 radians and continuous frequency in the interval (-f,, + f,), 
or even the interval (-~, ~). In theory the above model does not hold 
for the continuous case. The power spectrum for the infinite but count- 
able case would be a set of Dirac delta functions of amplitude a 
standing on the points (2, mp), f,) of a three-dimensional frequency 
space. The continuous case produces a power spectrum, S,(2,m,f), 
which is everywhere nonnegative and in general continuous over the 
region of three-dimensional frequency space where power is assumed to 
exist. A reasonable model for n(x,y,t) in the continuous case must be 
determined. Consider a single wave element 


a, Exp(i2m(2px + my + f,t + On)). (Qo) 
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The energy or mean square in this element is a2, Assume the ele- 
ment is a part of a egntinuum of elements for -~ < f < +~ and 0 < 6 
< 27m. In this case a, must be an infinitesimal energy associated with 
the frequency differential, df, and space frequency differentials, dt, 
and dm, which are related to the direction 6 of the wave element as 
before. Let the power spectrum, S(2,m,f), be defined with units of 
amplitude squared and divided by unit spacial frequency, 2%, unit 
spacial frequency, m, and unit time frequency, f. The power spectrum 
is then a spectral density yeas at (2,m,f). In this case we must have 
the infinitesimal energy, a, defined by 


ay Ss S (i,m, 4) ad adm AS (2.15) 


The real valued, nonnegative function S(%,m,f) is a power (energy 
density) spectrum of the standard type in three-dimensional frequency 
space (2,m,f). Intuitively, we can write an infinitesimal wave ele- 
ment as 


[ext (é am (4x + mY vit t «)) VS(4,. 3.) AN dm 44 x (2516) 


where the positive square root is assumed. To arrive at a noe of 
n(x,y,t) for the continuous case, we need only form a triple "sum" of _ 
the infinitesimals or, to be precise, the triple integral N(x.y.t) = 


Jf (lr Gen army oe earl) S (Rm) dd dn af (69 


For different sets of values of the phase relation a(%,m,f), the 
wave system, n(x,y,t), has a different shape, even when S(2,m,f) is 
fixed. In fact there is a wide range of possible shapes of n(x,y,t) 
for a given S(2,m,f). The above development is more intuitive than 
mathematically rigorous. ey has been shown by Pierson (1955 - pp. 126- 
129) that if 27a(2,m,f) is a random function such that for fixed (2,m,f) 
phase values of the form 21a MOD 27 between O and 21, are equally 
probable and all phase values are independent, then Equation (2.17) 
represents an ensemble (collection of all probable) n(x,y,t) for a 
given S(2,m,f). The random process represented by the ensemble is then 
a stationary Gaussian process indexed by the three dimensions (x,y,t). 
Detail discussions of the above can be found in St. Dennis and Pierson 
(1953 - pp. 289-386) and Kinsman (1965 - pp. 368-386). The fact that 
a particular sea-way can be considered as a realization of a stationary 
three-dimensional Gaussian process has been verified. Refer to Pierson 
and Marks (1952). 


The model in Equation (2.17) as a stationary random process will be 
assumed in following discussion. 


3. A DIRECTIONAL WAVE SPECTRUM 


The power spectrum S(2,m,f) is the directional wave spectrum of 
n(x,y,t). If n(x,y,t) is to be real, every infinitesimal of the form 
given in the continuous case above presumes the existence of its com- 
plex conjugate. Let us consider the one-sided power spectral density, 
S'(L5,M5,£,)5 of a single real wave where 0 < f, < ~. For such a real 
wave element of length i,, from a direction 9), ORO <2 the value 
8" Qeotlontiolo S Si@veamapee nas S(-25,-m,»-fo)> where £, = K, cos 8; 

Mm, = K -gind,, and Kg = 1/A,. If the above real wave came from a 
direction (27-6), the power density would be Si Cessip te) = S(£,,-m,»£,) 
+ S(=253Mp»-fo) - 


Figure 2 illustrates a real wave of length A>, from a direction 6, 
in two-dimensional spacial frequency (wave number) space. 


Wave Direction 


FIGURE 2. REAL WAVE IN WAVE NUMBER SPACE 


If the wave number relation 


Ka 27 ¢* /yTanh (at Kh) (Saal) 


holds, refer to Kinsman (1965 - p. 157) and Munk et al (1963 - p. 527), 

where h = water depth, g = acceleration of gravity, and K = wave number 

= 1/\, being the wave length; then a relationship between f and (%,m) 

is implied that requires a wave frequency fy to have a unique wave num- 

ber k,. From this we have the general one-sided spectral form for waves 
where f = f, of 


2 


2 


i] 
ta 
N 


S' (2,m,f£,) = (zero where g2 + m2 # Ky 

a power density > 0 for g2 + M on 
S'(2,m,f,) thus defines power density at f = f, for wave energy over 
0 <6 < 2n. Figure 3 illustrates this case in wave number space. 


S(L, mM; 385) 


Sis (lapap may see) 


FIGURE 3. DIRECTIONAL WAVE SPECTRUM AT A FIXED FREQUENCY, fy 


We want to estimate the shape of S'(2,m,f,) above the circle 22 + m2 = 
K_ “ in a directional wave train analysis. Remember, the S'(2,m,£.) 
above is restricted to fo > 0 and is, in fact, equal to 


t S(Xm4) + S (-R-m, -4.)]) 
S (2h, m4.) = §(- A,-m,-4,) 


where 


(3.2) 


if n(x,y,t) is to be real. 


Let us see how S' (2,m,f) might be found: we have said that 
n(x,y,t) can be assumed to be a stationary Gaussian process. One char-_ 
acteristic of such a process is that for fixed values (x9, yo, t=) ota 

_the process indices, n(x9. Yo, to) is random variable with a Gaussian 
distribution; i.e., 


Ne -p\* 
Preb (n (%0. 4, te) < n)= Sas gt (F*) dy | 


where yp is the arithmetic mean of n and o@ is the variance. Intui- 
tively n is as likely to be positive as negative, so let us assume that 
Prob (n(X9, Yoo to) < 0) = 1/2. Since n is Gaussian distributed, and 
is thus symmetric about its mean, we have Prob(n(x9, yo, to) < u) = 1/2 
or that » = 0. For o%, we have (using expected value notation) 


gta El(r-a)] = El (4.2) 


where we are thinking of n as a random variable. 


A Gaussian process is completely defined statistically if we know 
the form of the mean 


E(n(x,y,t)) and the covariances 


- [x (Kat) - EO (x4, €))] q(x +X, +, t+T)- 
E (» (xeR yeY,t +T))]} : 


where X, Y, and T are space and time separations, respectively. Refer 
to Parzen (1962 - pages 88-89). We have assumed E(n(x,y,t)) = uy = 0 and 
that the process is stationary (only weakly stationary is necessary). 
Hence, by definition of weak stationarity, we have for each (x,y,t), 

and yet independent of the particular x,y,t values, the covariance form 


RYT) SE Gst) Qkex yey, teT)] OD 


All of the properties of the stationary Gaussian process n(x,y,t) are 
implicit in R(X,Y,T), just as a knowledge of uy and o“ for a single 
Gaussian random variable completely defines such a random variable. 
Here it is important to understand that we are discussing expected 
values across all possible realizations at a point (x,y,t); i.e., 
across the ensemble of all possible sea wave shapes at (x,y,t) for a 
given S(2,m,f). 


There is a simple and unique relationship between R(x,y,t) and 


S(2,m,f£). Consider a single real wave element (from Equation (2.16) and 
(3.2)) as a random process and write n(x,y,t) = [EXP(i2n7(2x+my+ft+a) ) 


+Exp(-ian(Lasmy +ft + a))]Ve my AT dm df. 


Form the covariance function 


R(XY T) = E}(¥.¥.4) n(s+X.4+Y, t+T)| 

= Elexp (cen (9 (2x+X)+m(2u+Y)+4(et+T) +20) 
+ Exp (2 2a (M(axeX)am (24+Y) +f @t+T)+20) 
+ Exp(iat ({X+ mY +fT)) 

+ Exp (ian (AX+ mY¥+ #T))] S(t.mf) dt don Af 


where E is the expected value over the ensemble for any fixed x,y,t. 


Note that 
\ S (im. f) dhdm df 


is a constant with respect to the expected value. 


Consider the following problem. Let u be a random variable with 
uniform probability density function 


Fe K OS US 2 radians 
fe) else where 


Define a random variable Z = ell, 


as 


The expected value of Z is defined 


(em f(uydu 


ott 
K(edus 
6 


Ete) 


Considering the random variable nature of the phase 21a as described 
following Equation (2.17) at the end of Section 2, and applying the above 
notion to the cross product terms of R(X,Y,T) we obtain 


R(XYT) = [ exp (izw(lXe+emy + *T)) + 
exp (-i2n(I1X+ mY +£T))]S(hmf)ad am de. 
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For a real wave element we have, where S(2,m,f) = S(-%,-m,-f), see 


Equation (3.2), R(X,Y,T) “Exp (é 2n (IK mY+ £T)S (Son, f) dl dm df 


+ Exb (-izw(QX+ mY + fT )) § (-4-m,- f) ALaAmdf 


which is simply the sum of the covariance functions of two complex 
wave elements which are conjugate pairs. It also follows that 
R(X,Y,T) is real valued. 


Reverting to the complex wave element form, and noting that the 
expected ensemble value of cross products between different wave ele- 
ments is zero in a manner similar to the case of cross products shown 
above, we obtain the composite general relationship 


R(XYT\= f(( (cep( 20 (IK + mY + {T)S(dm.f)dtdads C4 


eo -63 ~ GD 
We have demonstrated, but not rigorously proven, that the covari- 
ance function R(X,Y,T) is the three-dimensional Fourier transform of 
the directional power spectrum S(2,m,f). 


We cannot hope to be able to estimate R(X,Y,T) for continuous 
values of X, Y, and T. However, there is a way around this problen, 


we can write the above as 
R(XN.T) = 


2 60 


fexp (iawfT) (exp (iate (2X mY) S(im.f)Mddeld$ 5.5) 


i-@® = 60 


which is in the form of a single dimension (variable f) Fourier trans- 
form of the term in brackets [ ]'s. Note this term is not a function 
of T. It depends only on the value of (X,Y). Further, by Fourier 
transform pairs we can write this expression as 


[ Ie = ( (xysr) exp ten fr) a7 (3.6) 


In general, assuming that the term in [ ]'s is complex, we can write 


[C(KY,F) -2 Q (XN F)J = 
J R(xyr) Exp(-a2ufT)AT | (3.7) 
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To find [C(X,Y,f) -— iQ(X,Y,f£)] we need only know R(X,Y,T) for continuous 
T for the given value of X,Y). Further, we have just stated that 


cc(ays)-ioayf)l = [ ].= 
(Sit.m#) Exp(i2t (AX+mY)) dam f (3.8) 


This is in the form of a Fourier transform; thus, we can write from 
transform pairs 


Stns) = ff[COKY,4)-0 0 (KY F) exp (Lem (1K mY) Axa 


(3.9) 


As has been stated, we cannot hope to have a continuous set of values 

of (X,Y). The solution is to find R(X,Y,T) for continuous T and 
selected values of (X,Y), and then employ the above to estimate S(%,m,f). 
This is described in the next section. 


4. CROSS SPECTRAL MATRIX OF AN ARRAY 


Let us look at n(x,y,t) at two fixed points in space, say 
(X9 Yo) and (x ,, y,). This would give two stationary Gaussian pro- 
cesses indexed on time alone because (x), yg) and (xj, yj ) are fixed. 
Thus, we may write 


Holt) = 4 (Xo,40,¢) and 4a (t)= 4 (x,y,6) 


If X = (x, - x9) and Y = (yj - yo) Then we can say, since n(x,y,t) is 
assumed weakly stationary (see Equation (3.3)), that 


R(XY,T) = Elv(t)-: % (t+T)] 


where the expected value is over the ensemble for some specific value 

of t, where - ~ < t < ~, Let us extend this idea by a change of notation 
and let N(X,Y,t) be a two-dimensional (vector) process, double-indexed 

on time; i.e., let N be a vector function 


N(XX,t) = (7eft)s 4 (8) = N (€) 


for —ce<t< oa oe) 
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We can then write a generalized covariance function (assumed to 
be finite) as the matrix equation 


R(T) = E[ Ne) N(t+T)| 


= JE (no(ed acer t)) E(re(#) U (t+ 7) 
E (a(t) nle+T)) Eln,@) le +T)) 


R (0,0,T) R(XW,T) 

T | Rtx-y¥7) R(0,05T ) (4.2) 

Now R(-T) is 
R (0,0, -T) R(x,Y,-T) | 
R(-T) = 
a R (-XNX-T) R(o,0,-T) 

and R(-T) transpose is 

s Rikon.) Rien 

RCT) 

R(X, Y,-T) Rio, 0, =F (4.3) 


We have by Equation (3.3) and stationarity that 
R (-X,-Y¥,-T) = E (% (A, 4 t+T) 4 (Ke, 4e,)) 
= E (7 (x0, 40, #) 7 [x,,4,,¢+T)) = R (X,Y,T) 


(4.4) 


It then follows that 


CG ae R'(-T) (4.5) 
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From Equation (3.7), 


ee PRY. #) = [CONF)- LON] 


and we get ie 
PUXY,£) = [ RIXY,T) Exp(-Lews TAT 


Now R(-X,-Y,T) = R(X,Y,-T) by Equation (4.5). Thus, we have from 
Equation (3.4) 


R (X7,-T) = {((expli an (AX + mY + (-T))S(DomF)dldondt 


or, by the same procedure, that Equation (3.6) was obtained from 
Equation (3.4), we have 


P*KY,#) = [ R(X,Y,-Tlexp(i2w4T) AT (4.6) 


Now, since R(X,Y,-T) is real and Fourier transform pairs are unique, 
we must have from Equation (4.6) that 


fr (X,Y, -T) Exp (-i2t £T) 4T =P (X,Y, #) 


a parallel form of Equation (3.6). Therefore, the Fourier transform of 
R(-X,-Y,T) = R(X,Y,-T) is the complex conjugate of the transform of 
R(X,Y,T); i.e., (see Equation (3.7)) 


P(-x,-v,f) =[PIXYO] = PY 


In general, we have that the Fourier transform of R(T) = R! (-1) is 


Ploo,f) PUXY.F) 
P(X,Y.4) P(o0,0,¢ ) (a7) 


Note: P*(0,0,f) = P(0,0,f) is real. 
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By definition it follows that 


CEGYG 2) = RelPi(xy5£)1 


Q(X,Y,f) Im[P(X,Y,f) ] 


where -~ < f <~ , 


The function C(X,Y,f) is called the cospectrum and Q(X,Y,f) is called 
the quadrature spectrum. Both are spectral density functions. 


Explicitly, we can think of (x, y,) and (x , y ,) as being the 
location of elements of a probe-array with space separation (X,Y). 


The first step to find (or estimate) S(2,m,f) (see Equation (3.9)) 
is to find P(X,Y,f) related to a pair of array elements. This is a 
problem of estimating the cospectrum and quadrature spectrum of a two- 
dimensional (vector) stationary Gaussian process. Goodman (Mar 1967 - 
Chapter 3) has an excellent treatment of this subject, which we will 
discuss. Kinsman (1965 - Chapters 7-9) also discusses the subject. 
The essence of the problem is that if P(X,Y,f) is continuous and negli- 
gible for | £ | > f£,, then R(X,Y,T) can be obtained by a time average over 
a particular realization N(X,Y,t) instead of having to average (find the 
expected value) over the ensemble. This says that we can find R(X,Y,T) 
by obtaining two time series (realizations) r,(t) and r;(t) measured 
over time at only two points; e.g., (x9, yo) and (x1, y1) where 
(x] - Xo) = X and (yj - yo) = Y. The relationship between (x, (t) r, (t)) 
and R(X,Y,T) -~ < T < ~ is 


ti 
RIX.Y. T= Limit & | wlt)en (t+T) dt 
() ty 


(4.8) 


where ro(t) is a realization measured over time at (x , y_) and r4(t) 
is measured at (x1, y,). We can simplify the notation by an expression 
for a time average given by 


ROG) = RG CL) = PACE) GEED) i 


It follows from Equation (3.7) that C(X,Y,f) = Co,(f) = 


po 


| Ro1 (T) cos(2nfT) dt, 


-co 
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and 


ios) 


Q(X,Y,f) = Qo1(£) = il Ro j(T) Sin(2n£T) dT. (4.9) 


-& 


We also have (see Equations (4.2) and (4.7)) Co (£) = Cy,(f) = Po (f) = 
Paice) 


Qo0 (£) 


P§, (f) S Coy (£) -i Qo1¢£)- (4.10) 


Qin Ce) 40); 


The phase of Po (£) is given by 


SOLES (4.11) 
Co Cf) 


Qq3 (£) = Arctan 


This is the expected phase lead of the signal at (x 
sienal at (Gol ya) for £ where) —)o<) f o< 27. 


°°? We) over the 


For an array of N depectore located at (xj, yj), Cx y2); (x3, Y3 ) 


» (xy Yu) we can find N eee te ] 202, N,, die N. This nee 
a Srone spectrum Pai (Pdi S oo and since Pij = 3, (see 
Equation (4.6) and 4.7)) 
we have 
N@=) (4.12) 
2 


unique cross spectra or a total of [N(N-1) + 2]/2 unique spectra. Thus, 
for real n(x,y,t) we have Pij (£) = P¥,(£) and Py (£) = P# (-£) allowing 
us to define the information” about Pry, ¥0538)) obtainable alln an array by 
a cross spectral matrix 


mm) Gr) eo = = Corl) 
QU) P.(4) oe Can (Sf) wnere 0¢ § < oo 


Q..(4) eae Pua) (4.13) 
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This information can be used together with Equation (3.9) to obtain an 
approximation to S(2,m,f). Numerical details for finding the spectral 
matrix are given in Bennett, et al (June 1964). 


It should be pointed out that negative frequencies are still con- 
sidered in the relationships being discussed. We do not know P*(X,Y,f) 
for continuous values of (X,Y). We do know from the spectral matrix 
the values of 


PT yf) fer Cates be begn 


where X.. Ge x) 
at j i 


(y. - y.) 


ij j aL 


We also have from Equation (4.7) that 


VP Cig tod) = POs ip Yi @a8) « (4.14) 


From Equation (3.9) we get 


S(tim.f) = ieee Y, ) cos(2m(&X+my) 


~iPTX.Y, 4) sin [211 (K+ mY] PAX AY 


or, since S(£,m,f) is real 


S(t.m.¢) = (( te (xy £) cos[e2T(Ax+mY)] 


-ao —~oO 


-Q(x Y.)sin[en(1X+mY)]} dxdy 


(4.15) 


Let us consider treating the points of (X,Y) where we know C(X,Y,f) and 
Q(X,Y,£) as weighted Dirac delta functions; e.g., at (Ky o> Y42) we get 


C (X,Y, 4) = be KarXasf) F (X- Kn) df (Y- Yel - 


IL7/ 


Reverting to the Ci, Qi5 notation of Equation (4.13), we have, where 
Xi5 = -Xj; and Yj; = -Yij, C54 = Cis and Qj4 = -Qij- The numerical 
form of Equation (4.15) then becomes 


-| 
S(£,m,f£) = by, Cyy(£) + a= S bs} Gy cos [en (tx, + mY.) 
=! J=itl 
Choice of b;, values is arbitrary. A reasonable choice is b.. = 
[N(N-1) + 1)71 for all (i,j) (refer to following section). We now have 


a basis for an approximation of S(%,m,f). Before exploiting this result, 
we need a few side results. 


5. SPECIAL CROSS SPECTRAL MATRICES 


Assume that we have a real sea wave of frequency f_ > 0 moving 


from direction 6, where 2, = K,Cos6, and m. = K,Sin 85, Ky being the 
wave number from Equation (3.1). We can write the wave as 
n(2%,4,t) = Acos(2w(Lx+mey +f.t +a)) (5.1) 


Since the root-mean-square (rms) value of a cosine wave is A/v2, we have 


for the two-sided (-»~ < f < ~) directional power spectrum of the wave in 
Equation (5.1) 


S(Xm,f) = A [S(1-2) Sn-m.) SEA) + 
+ §(%40.) §(m+m.) S(f al (5.2) 
or in polar form where K = /y2 emer Gi= arctan (7) 


S(k.e.4)= A [S(e-0.) S(#-£.) + 


+ §(0-(e.-m))d (+f) ] §(K-K.) 
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From Equations (3.8) and (5.2) we have for a single wave of frequency 
f, > 0 from a direction 6, that the two-sided 


Ret) = A [exe (ien(t 0X, m, Y, NS (F=f) + 
+ exe(iat(-d Kom ¥;))4 (4 + f.) | 

or wHere ne) = C4) — £Q, (8) that 
C(t) = Cil-f,) = ne 
C= C6) = B cos (2m(t, X,0 m.¥,)), 


Q(t) = - Gy (-4) = — 4 swn(en v( Xm ¥) 6.9 


describe the elements for the spectral matrix of a single wave. Recall 
that, in general, Ci; (£) = C54 (£) and ay (Ge) SO) Ge Substituting 
into Equation (4.16), we get where bij [N (N- DET 1/M (5.4) 


AX 
and f, > 0 is assumed S ) m, £ At A 
: U mf) = 


Kl A 
2 cos (20 (1. X,+ m%;)eos (21 (LX,+m ¥, )) + 


+sin(2m(1.X,+ me¥;)) SIN (271( 2X; mY,))] oR 
S(bmf.)= 


Mh 02d, 5 cos (2n[ X(t ~{.)+¥, 3 (m -m)] (5.5) 


seiel 


N(N 1) 
2 


>’ 


w=! a 
The choice of bij =ely/ Mi [N(N-1)+1]72 and observing mick & 4 
Seo Jsley 
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gives the convenient result (Note: 
in place of (25M) ) 


S) (L., ma, f ) = S(- JN. ,-me ,— fe) 


=A{i- Sen = 4 


FORps toms 0 we would use (ages) 


4m 


which agrees with the values of Equation (5.2) at the points (%,,m ,f)) 
and (Boosey Note that there is not general agreement elsewhere. 
In fact, Equation (4.16) may (and does) give negative values for the 
approximation of S(2,m,f). This is a problem of probe array design and 
is directly related to the directional resolving power of a probe array. 
This problem is discussed later in another section. 


Consider the case of a single frequency, f, > 0, real sea with 


equal wave energy from all directions; i.e., isotropic waves. In this 
case, assuming the wave equation (Equation (3.1)) holds, we have 


S(t,m,f)= [A §({0'+ m=) - (97 + m2)) ][s(¢-t.) + 4 (+4)] 
Savonene Kas Q"+ m?, K=+\k* , ano 6 = aRcran ("%) 


S(K,0,f)= AS (KEK?) [i (F-f.) +5 (F442) J 


(5.6) 


as the two-sided directional power spectrum for such an isotropic wave. 


Since % = Kcos® and m = KsinO, Equation (3.8) can be expressed in 
polar coordinate form as 


WV oo 


* 
Bids | IS(K.6,4) exe (i 2m K (X,cos 8 +Y,5IN 8) KAKA® 5. 7) 


“nu 


Dace 2 
Letting Djj= (Xi + XG 
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and g.. = ARCTAN (| 
4j Xi 


we can write 


aye W co 
BUF) = § [S(Kye,#) exe (Lem kd,jcos@-Q))Kdkda 5 
-—W Oo 


Using Equation (5.6) for S(K,9,f) we get 


Bld =A K [exe (Lam K,D,cos (6 - o,.)) de 
Lat 


tf 
= (ih K. [cos (at k,D,jcos (8 —$,\\ de + 


=“ 


+iA K. [sin (2% K, D,,cos (@ -4,))d¢ 


-1 


(So) 


Now, departing from the above development, consider the following 
integral where 


z= ew Kk Di, = 9; 
AT 


(cos (n9) Cos (z cos(e-w))de 
-w 


Let > = 0 — py and we get 


(cos ing +nw) cos(z cos) do 


Zi 


Expanding cos(né + ni) we get 
Wy 
cos nW cos np cos (zc0s ) do 
-W- 
= 
— SIN mip {sured cos(zcosh)dd . 
n= 


We have (rt - wv) - (-n -¥) = 27 so that the integrands are over 27 
allowing us to write the above as 


mig 


cos nw cos ng cos(zcos 6) do 
a 


aT 
—sinnp) sin ng cos(zcos >) do . 
= 


From Ryzhik and Gradshteyn (1965 - page 402) we have (since sine is odd 
and cosine is even) the result 


(cos (ne) Cos(z cos (e- wp))d ") 


= cosnp|an cos (22\ Jz] , (5.10) 


In a similar way we get 


T 
{cos (n6) SIN(2 Cos (O- p)de 


<1 


= cos np[2t Sin es) Je] : (5.11) 
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where J, (4) is a Bessel function of the first kind of integer order n. 
Returning to Equation (5.9), the above results give, for n = 0 where 
ee 0, the result 


R (tt) = Aen, Je (ank.D,)  - 


Since the above is real 


C,(#$) = 2mkA 
Cet) = awk AJ,(em kD.) 
Q,(¢ 4) == 0 


(5.12) 


describe the elements for the spectral matrix of a single frequency 
isotropic sea. 


The above two special cases for sea waves are the extremes of 


directionality of real sea waves of frequency f, > 0. These results 
will have important applications later. 


6. A MEASURE OF ARRAY DIRECTIONAL RESOLVING POWER 


From Equation (3.9) we have the Fourier transform 
S (tm. f\=((P'(x,y, eexe(—Lem(Ux+mY)\dxdy «1 


In practice we use a probe array with elements at (xj, y;)..-, (X> yx) 
to obtain P*(X,Y,f) at the separation points (0.0), (Xo, Yy9)> 


(-X19 »-Y49) 2.900 (X-1,k> Yea ke)? (-Xx-1,k> Ge aL ie) 


We do not then know P*(X,Y,f) but the product 


PEC) CGN) (Gm) 
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where g(X,Y) is a set of Dirac delta functions standing on the separa- 
tion points of the probe array, and zero elsewhere. 


Thus, we have the estimate 


S (wi, f= (Pr, Y,4\ at Y)Jexr(- Lan(ix+ my))dudy 


G (t,m)=f (gx) exe (-ien(4X+mY)) dx dy 


-@ = 


(6.3) 


Using this and Equation (3.8) we have for a given (%,,m),f,) that 


GO, 62g, & 


8 (Lm, f) = f \ (Sm) ene (iem(Ik+mY)) dtd] (ay) 


exe (-i2m(tX+m¥))dXd 


= {f(s S(L,m,f.)9(X,¥) exe(-aetr(x (1. ~\)+¥ (mo- -»\)) xd} 


=60 -@ ii 


={(S(0. m. offal 9(X,YexP(-ien(x(2.- -1)+ (m,-m))) axa yfadn 


= [(S(tm.4)6( Vale m,-m)dddm 
Brees (6.4) 


As expected, $(2,m,£) is a two-dimensional convolution of the true 
directional spectrum S(2,m,f) with G(2,m) the Fourier transform of 
g(X,Y). We see then that S(2,m,f) is a weighted average of S(,m,f) 
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and that G(%’,m) is a measure of the directional resolving power of the 
assumed probe array. By the nature of g(X,Y) we have from Equation 
(6.3) 


G (im) = 1428, d cos(2n/ (Lx;+ mY,;)) 


(6.5) 
! 


If we assume that the wave eqyatiog (Equation (3.1)) holds, we find 
that S(2,m,f,) is zero when 4” + m* # ae the wave number for f, 


From this we have for a given (29M) that the directional eRilelae 
power, DRP, is 


DRP (A, mf.) = (a(n mm) — Ks) G(I.-0,m.-m) Addm 


-a@ —@ 


2 > K 
rer h=k,.cos@ » M= Ke SIN O, WHICH IMPLIES Vem =K, 


and we get, for energy coming from a direction 6, at frequency f, as a 
function of 0 < 6 < 2n, that 


DRP(@ | @.,£.) = G(k.(cos @,—cos@), K.(sin@,—SIN Q)). 


From Equation (6.5) we get 


DRP (e|¢.,f\= = 1425 deos[en( X,jKe (cos@,—cos €)+Y, |Ke 


d=) grat! (SiN 0, — Sit é))| 
= 1425 Scosfen (Ch -k,.cos 6) X, + (n,— Ke Sin 6) Y, “I 


4=i JzAti 
Nore THAT DRP(0,| 0, £) = N(N-1) +1 
Compare Equation (6.6) and (5.5). Except for the amplitude tern, A2/4M, 


in Equation (5.5) the equations give identical results. The choice of 
= [N(N-1) + 1]7~- = 1/M is again found convenient. 


(6.6) 
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7. DIRECTIONAL ANALYSIS FROM THE CROSS SPECTRAL MATRIX 


First, we consider a fundamental approach. We have for a pair of ; 
detectors I and J located at (x;,y;) and (x,,y,) respectively, the 


cross spectral matrix, P .(f) = Cij (£) + iQ;.(£) and more impor- 
tantly ¢14j(f£) the phase thad @ie I ae J etal by, 


O(f) = ARCTAN foun 


C;; (€) ‘ 1) 


The actual phase lead may differ from this value since the true phase 
lead 68 is some one of the values 


Q(t) +hat 


where 
Ri=0,£1,42, 


Consider a single wave of frequency f, with corresponding wave 
length \, and wave number K, = 1/\,, and find the direction the wave 
must travel; ioGog alia 2 etaele wave to the spectral matrix results for 
the detectors I and J. 


Let D, ij be the distance between I and J. The distance between 
I and J in wave lengths is K oDij- In radians this is 2nK.Dij;. From 
this relation we get 


—2Tf K. Dy BS d< eT K. Di 
or that only values of h such that 


—2uk.D (A +hew S ett K, Dij G2) 


give physically acceptable candidates for the value of >. If Dij < A/2 
only one h value is valid. If Nall 2 < Dij < ro at most two values of h 
are valid, etc. 


The problem is how to find the true direction, 9,5, of a single 


wave given the above possible value(s) of ¢. Consider a given value 
of » in terms of wave length units and we obtain 
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ink ne “ork, * (7.3) 


Figure 4 illustrates a case for > 0 (and thus Lg > 0). From the 
figure we have, where the true direction is 84 and true phase is 9, 
that 


6, = Thank ats +H 


nol = (7.4) 
WHERE S\ = ae KD i 
oR “a= arcein|—_O _ = arcsin| D | 
ay K, Dij kj 


FIGURE 4. WAVE DIRECTION ANALYSIS 
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Recall that sin(a) = sin(+7-a). Thus for a given > > 0 we get, since a 
must be obtained from an arcsin relationship, two possible values of 
wave direction, the true value 6, and its image 


es ee at bad 
= yj- tt a= OG) == Wii- Loe a | (aD) 


This is illustrated in Figure 5. In actual practice we do not know the 
true direction, thus a given value of > > 0 gives the direction as 


a= Wi* ie ec] 


where 0 < a < 1/2 is obtained from the principle value of the arcsin. 


If ¢ < 0, then I actually lags J by | ¢ | > 0 and Lg = Lo/Dij <O, S© 
that arcsin (Lg) gives -1/2 <a < 0. 


True Wave 
Direction 


8 


° 


wy, Q@’ (False Direction) 


FIGURE 5. DIRECTION ANALYSIS FROM A PAIR OF ARRAY ELEMENTS 
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The implied directions, for » < 0, are directly opposite from those 
for ¢ > 0 so that for a given value of 9» < 0 we get directions (a < 0) 


or 0,= at = + Ow 
at 
a 


as before. 


Thus, where the principle value of arcsin is assumed, we get a set of 
possible directions 


8,, = Wi elon + 94] 
h, = ARCSIN Gil)+ her 


aw K, Daj (7.6) 


where 


h being constrained by 


Qii( ‘eben an 
An example of this type analysis, for an array pair, is illustrated 
in Figure 6. Thus several estimates of 6) are available (at least two). 


The estimates of true direction, 6,, often vary from one array 
element pair to another, making the selection of a true 6, value diffi- 
cult. The selection is also hindered because half of the estimates of 
84 are of the image type; i.e., false estimates. 


While the above directional method leaves something to be desired, 
it does illustrate the basic directional information produced by an 
array of detectors. 


A better method, suggested in Munk et al (April 1963), of using the 
spectral matrix directional information to fit a single wave at each 
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FIGURE 6. DIRECTIONAL ESTIMATES FOR A PAIR OF ARRAY ELEMENTS 


frequency is given below. It is based on Equation (4.16) in the 
form 


S(tmf) = £ (CoC) + ao y L[ealtdeosten (Lx4j+ mY) 


(a7) 


cae +m Y,, mI 


where M= [N(N-1) +1], C(f)= = 2 C..(f) 
=| aA 


and N = number of array elements. 
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Recall that for a single real wave of frequency £, > 0 and known 


direction 6, Cy 4(f ), Q;.(£,) are known (see Equation’. (5.3)). These 
values give j 


a 


S{k,.6,.4) = S(k,,f, -W-aiy=2 


When a single, well-directed swell is expected, it is reasonable to 
assume a single wave for a given frequency exists, and to select 6, and 
Ao = A2/4 such that the least square error between the Haeeraieleaill, 
cross spectral matrix for a single wave (see Equations (5.1) and (5.3)) 
and an observed cross spectral matrix is a minimum. Accordingly, using 
the expressions of Equation (5.3) and an observed cross spectral matrix 
for a given frequency, f,, we can form the squared error 


= (C, -'A,)? 


i 22 [6;- —A,cos(at(h. X;,+ mN;, MN 


+203 (0, a Assin (an  hXpce mY a 
ist jet 


isl it (7.8) 


Expanding and Matte terms we get 


H= Sone (Cis +Qi3) + [NIN=1) +1] AS 


431 Joie 


~ ads eye Cis cos(am(te X,,+ m.¥,,)} 


is year 


nO jo!N (am ( (4.%; yt Me x) (7.9) 


AZl Jeaay 


or using results in Equation (7.7) 


Syil 


H= tea Dl Ci +Qi5) +[n( n-i)ei](R-2A.Sllomt] 


Azi j= Lt 


To tind A, that minimizes H, consider 


= = 2(n(n-i\+ iI[A- = SID), ‘| =0 (7.11) 


° 


This requires that A, be of the form A, = S(2,m,£,) and a resulting 
value of H of the form 


=) oD) (C +95) ) = [n(n -1) +i)” 


is) jties (iat) 


Since Co. Cra and aan are all nonnegative, a minimum H results when 


Si@vemts ES) is a maximum. A choice of 2, and m, that maximizes S(2,m,f£,) 
mreiecleens a 0, = arctan my [Xo which is onieinemn. Remember that we are 
assuming Ko = y2 + Me holds. along with the wave equation. The results 
then for aeen fo is a two-sided energy spectrum estimate A (£28 =) 
Appendix B ORIEHINS a listing of a FORTRAN II program for finding 

ACES so) > the least square wave fit, from a set of spectral matrices 
obtained from the task SWOC data collection and analysis system 
described in Bennett et al (June 1964). 


Examples of least square single-wave train fit analysis from 
Bennett (March 1968) are shown in Figure 7. 


A more complete collection of the directional spectra calculated 
from data collected off Panama City, Florida, is given in Appendix A, 
and in Bennett (November 1967), and Bennett and Austin (September 1968), 
an unpublished Laboratory Technical Note TN160. 


We would actually like a continuous estimate of S'(2,m,f). Con- 
sider then a third method. From Equation (3.8) we have for a pair of 
detectors, as illustrated in Figure 8, that 


(Yt) = [(s (L.m.f) exe [ian ({k+mY)] didmn 


-2 —@ 
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FIGURE 7. LEAST SQUARE SINGLE WAVE FIT DIRECTIONAL SPECTRA Aj e))) 
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(Xp ’ Yo) 
Detector 2 


ro (X15 Yr) x 
Detector 1 


FIGURE 8. DETECTOR GEOMETRY 


Let 2 = K cos@ and m = K siné 


D= ix? at: ye wv = Arctan : 


and 
or 


X = D cosy and Y = D sin bp. 


Using the above changes of variable, we get 


M co 
P"™(x,y,¢) = | [stk af) EXPlien (K D cos 8605 H+ KOSINOSIN 9) 
-% “Oo Kdkde 


(7.13) 


Prixw.s\=| (3 (K0,f) Ex pliant KD cos (4- y)kaK da 
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We have from Equation (3.2). that 
S(K,e,4) = S(K, 6-1, -f) 
If we think in terms of fp > 0 
a 
Ss (K,6,,) =e S (k,8,£) a) 
We are assuming that the wave number relation of Equation (3.1) holds. 


Thus Figure 3 is applicable, and we can write the one-sided spectral 
density as 


$'(k,e,f) = 2 (8,8) S (k-K) 
ich ( K= kK, 
wHerE (K-K,) hes Baie. 


This allows us to write, mao<¢ F< +eo , 


Rye te = { {a (6, f.)d(K-K,)exPliewk Dcos(o- WKdkde 
- O 


a 
PUXN fe (a (o,f) expiant K,D cos(o- yk, de 
-T (7a) 
We have reduced the preblem to finding [a(0,f, - K,]. 


From Equation (7.15) we see that 
P*/o,0,f)= P(f)= falet)Kde 
sd 


where P(f,.) is the power spectral density of frequency £,- For better 
comparison of cases where P(f;) = P(£,), [£,| # Te, | it is convenient 
to express a(@,£,) in a normalized form 


A(6,f,) = a(8,f,)K, 
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where we get 


a 
P(£) = JAle,f)de (7.16) 
-T 
Thus if the energy distribution as a function of direction is the same 
for P(f1) = P(f,) then we also get 
A(@,£1) = A(6,f,)- 


Consider now (assuming A(6,f) can be so expressed) a Fourier series 
expansion of A(6,f) for fixed f. Clearly it is periodic in 6 with 
period 27. Thus for any given f = f, we can write A(6,f) as 


A(9) = = +2[0,,cos n@ +b sin n6] (7.17) 


Substituting this expansion into Equation (7.15) we get 
S Vv 
* : 
ie (x,Y,*) = $ {xe Lav KD cos(6-y) de + > [a, feos ne 
24 
—-% w cs =F 
EXP LAN KD cos (6-y) de +b [50 ne EXP Lem KD cos(0-¥} 4] 
on P™ (X,Y, #) =f [ss 2m KD cos(e-) t)) de DAC Jasne cos(2tekD 
“Tt 


cos (e- “W)de +b of ne cos(2t KD cos(e-y)}e 
+ ile (os (at kD cos (6- y}de + 


oy (2, (cosne sin (21 KD cos(e-¥)de > 


N=! -v 


an 


ba [sm ne SIN (24 KD cos(e-¥) 0) (7. 18) 


36 


Thus we can express P*(X,Y,f) as complex infinite series with unknown 
CoecEflctents ag), zis) aie and) Dai by, -.- and constants defined 
by the integrals (tet Y =! 27KD: ny= (05 252,000) Of thertorm 


a\s 


J tos n6 Cos (2 Cos(@ -¥)) de (7.19) 
T 

(sw no cos(2 Cos (8-y))de (7.20) 
AE 

. 

(cos NO SIN(2 cos (o-y)}de (7.21) 
“h 

and 

sin no Sin(2 cos(6-y)\de (7.22) 
=m 


Consider Equation (7.19) where » = (6-) and dd = d6, \ being a constant. 
We then have on changing variables 


-\ 
(ces(nd+ni) cos(z Cos g) Ag 
ee n- 
= cos n ¥ Jcos ng cos (2605 6) dg (7.23) 
=N- 
1-4 
— SINT y (sin n @ cos(z Cos g)dg 
=v 
Since the integrands in Equation (7.23) are both of period 2n and the 


interval [-7-, 1-W] is of length 27 we can write the equivalent of 
Equation (7.23) as 


cos n v (cos nq Cos(% Cos 4) dg 


— 
— SIN ny {sin ng cos (z cos g) ad (7.24) 
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From Ryzhik and Gradshteyn (1965 - page 402) and noting that the second 
integrand is odd we get Equation (7.23) equivalent to 


cos ny (on cos (2) Je (z) | (7.25) 


where J,(Z) is the Bessel function of the first kind. Employing a simi- 
lar procedure for Equations (7.20), (7.21), and (7.22) we get 


P(X.) = Sef an J, (21 Ko)} + 
> [ances ny an cos (2) J. (am kD) + ba SiINnY 2% 
cos (3) Sn (an Kr) 


+4 [o. Cosny 2m ae (2m kD) 


nel 


+ by, SINK Y ev sw(nt) J, (2% K 0) | (7.26) 
Now nn ath a.) n odd 
cos (11 } = { (-l)"% even 
and N (3) = n even 
S\ oy ot gate 


Thus we get P™(xY,#) = Qe fam J, (en «0)} 
+2 [en J, (21 «0) (- la Cos nw + b,, SIN ny) 


4+6,° 


+4 (ae (anko)(-!) ye (a, Cos ny+b,snny)] 


NzV3,8°° (7.27) 
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From a spectral matrix of the form in Equation (4.13), M = N(N-1) +1 
different equations can be set up using Equation (7.27). This allows 
us to get a system of equations for any m of the unknown coefficients 
AQ, 445 425 ++. 3 by, by, b3, ... while assuming the rest of the 
coefficients are negligible. We can then solve for the m desired 
coefficient values. This has not worked well in practice for two 
reasons. The inverse of the matrix of constants obtained is sparse 
and often ill-conditioned. Further, if the wave energy is from a nar- 
row beam width (30 degrees or less), the first 100 harmonics in the 
Fourier series expansion can be significant. There is perhaps a 

more efficient orthogonal set of functions than the sines and cosines 
of the standard Fourier expansion. Search for such an orthogonal set 
should prove fruitful. One might start with Walsh or Haar functions. 
See Hammond and Johnson (February 1960). 


8. SUMMARY 


It is believed that the least square method of using the informa- 
tion in a spectral matrix is the best method presently available. 
Examples of such analysis can be found in several of the papers in the 
bibliography. A collection of ocean-wave induced, bottom pressure 
directional spectra from these papers is given in Appendix A. 


An iterative extension of the least square method can be found in 
an excellent paper by Munk et al (April 1963). Some details of this 
method are given in Appendix C along with an example result and a 
FORTRAN program for the method. 


There is merit to using the coherency, 


R.- (+) Bs Msc : 
: CRS Pee 


to form the weights Dag in Equation (4.16). Ome idea being explored is 


L. = hs a Ings 
5A ae pee ae 
+2 ines R 
421 42at! 
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APPENDIX A 


A COLLECTION OF DIRECTIONAL OCEAN WAVE 
BOTTOM PRESSURE POWER SPECTRA 


This appendix is a collection of the results of a least square, 
directional, single-wave train analysis of the cross power spectral 
matrix resulting from the analysis of ocean bottom pressure data. The 
data were collected at Stages I and II offshore from Panama City, 
Florida, during 1965. The data collection system and the estimation of 
the cross power spectral matrix associated with a set of data are 
described in Bennett, et al (June 1964). Augmented pentagonal arrays, 
containing six pressure transducers each, were located seaward of each 
of the stages. Stage I is 11 miles offshore in approximately 103 feet 


of water, and Stage II is 2 miles offshore in approximately 63 feet of 
water. 


Certain parameter values are pertinent to the directional analyses 
presented: the number of data points in each pressure data set is 
N = 1800; the sampling rate is once per second, At = 1 second. On the 
cylindrical polar plots frequency is the radial variable and compass 
bearing the angular variable. The vertical axis is logjg of power 
spectral density in inches*-seconds of water pressure. The frequency 
axis range is 0 to 0.3 Hz in 0.05 Hz increments. This is illustrated 
in Figure 7 of the report. In each plot title, the date, time, and 
location (stage) is indicated. The value WD is wind direction in com- 
pass degrees, and WS is wind speed in knots. Appendix B gives a listing 
of the FORTRAN II computer program used to produce the plots. ; 
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APPENDIX B 


A FORTRAN II PROGRAM FOR SINGLE-WAVE TRAIN ANALYSIS 


The FORTRAN II listing of an IBM 704 program for the least square 
single-wave train analysis of a spectral matrix obtained from an array 
of ocean wave bottom pressure transducers is included. The listing is 
from the FORTRAN to ALGOL translator of the Burroughs B-5500 and is 
syntax free at the FORTRAN II level. The mathematics of the single- 
wave train analysis is described in the body of this report. The data 
collection system and calculation of the required cross spectral matrix 
is described in Bennett et al (June 1964). The plotter subroutines 
GPHPVW and PFB3D are also included. Bennett (March 1968) describes the 
plotting technique in PFB3D which was used to produce the plots in 
Appendix A. 


FORTRAN TO ALGOL TRANSLATOR 
PHASE 1 FORTRAN STATEMENTS 


C LISTING OF AN IBM 704 FORTRAN PROGRAM FOR DIRECTIONAL WAVE ANALYSIS. 
C 13350"? C BENNETT JULY 1967 SWOC 78°301°8210 
Cc REQUEST=0268 
Cc LEAST SQUARE SINGLE WAVE TRAIN FIT AFTER MUNK 
DIMENSION IDC41)29FQ€100)2WNC100)»PERIODC 100) *WAVLGHCE 100)» THMAX( 100 
X)*DEPATN(100)2AVEP(100)2P(100»696)» EMAX(C100)*TEM( 146) 
DIMENSION 06656) sPS16656)5E(100572),THETA(72) »DEC73) »THZRO(1496)>» 
XH(100)2H1(€100)2H2(100) »G000¢6) 
DIMENSION DATA ¢€1000) 
E2TEN=0,43429448 
TWOPT=6,281853 
RT0257,29578 
DTR=0.0174532925 
FIV20.087266465 
REWIND 3 
REWIND 9 
CALL PLOTS(DATA(1000)»1000) 
CALL PLOT(0.92=*30.09°3) 
CALL PLOT(C2.592.5293) 
DO 5 Ieis7e2 
ZIsC171)*5 


5 THETACI)SZI*0TR 


10 


DO 10 [=1%6 

DO 10 J=196 

DC I»J)=0.0 
PSI(I2J)=0.0 
0¢1»2)=100.0 
0¢153)=100,0 
0¢194)=100.0 
DC125)=100,0 
D(1s5)=100,0 
0¢223)=117,558 
0¢224)=190,212 
0(295)=199,212 
0¢296)=117.558 
0¢394)=117,558 
D¢395)=190,212 
0¢3952=190.212 
0¢425)=117.558 
0¢4945)2190,212 
0¢596)=117.558 
PSY IS TRIG ANGLE 
D IS DISTANCE 
PSI¢€122)=0.0 
PSI(123)2=72.0#DTR 


PST(124)=144.0*DT 


R 


15 


20 


743 


22 


PS1€125)=216.0*DTR 
PS1(126)=286.0*DTR 
PS1€223)=126.9*DTR 
PS1€224)=162-0*DTR 
PSI€2#5)=198.0*DTR 


PST(226)=234.0*%DTR 


PST€3294)=198.90%DTR © 


PS1(325)=234.9*%DTR 
PS1(3%6)=270.*DTR 

PS1(425)=270.0#DTR 
PS1(496)=306.0*DTR 


PS1(6526)=342.9%DTR 


GONDC]1)= 1 FOR USABLE CHANNEL DATA AND 0 FOR BAD CHANNEL 


READ 2» ISKIP s(GDODCI)s» 12126) 


FORMAT(I2 »6F1.0 ) 
IFCISKIPI100230220 


DO 25 J=isISKIP 


READ TAPE 35TDsFQC 1) 2 WNC 1) 5MoK 


IF(K)24222021 


PAUSE 20202 


TAPE OUT OF PHASE WITH DATA READ DESIRED 


GO T9 15 
MP=M+1 


00 23 L=2»MP 


23 READ TAPE 3 
25 CONTINUE 

GU TO 15 

FQC1)=0.0 CAN NOT BE MEANINGLYFULLY PROCESSED 
30 READ TAPE 35 1DsFQC 1)» WNC 1) o4eK 

TF CK)21031921 
31 MP=M+1 

M0 32 L=2sMP 
32 READ TAPE 3,1D,FQCL)»WNCL) sMoK sDELTAT »DEPTH»PERIOD(L) ,WAVLGH(L)» 

XDEPATNCL) sAVEPCLI 9 COPCLelsJ)sJH196)21=196) 

P(L»sI#J)="0.0 IF CHANNEL [| OR J IS NO GOOD 

WNCL) IS WAVE NUMBER=2PI/ZWAVE LENGTH IN FEET 

PI=3,141592/ 000 

VALUE K NO LONGER NEEDED 

GOONCH=0.0 

DO 49 1=196 

GOODC1)=GOODCT)#PCAsToT) 

IF CGOODCI))41 240544 
41 GOODCH=GOODCH+1.9 

GOOD(1)=1-0 
40 CONTINUE 

TERMS=1,0+(GOODCH*(GO0DCH=1.0)) 

DO 70 L=2»sMP 


SUM=0-0 


DU 59 J=196 
IF CGOUDCI))51*50+51 
51 SUM=SUM4#PC(LsTrl) 
59 CONTINUE 
AVEP(L)=SUM/GOODCH 
DO 69 k=1»%72 
SUM=0-0 
D0 90 I= 195 
IFC(GOUDCI) 2909094 
91 1P = I44 
DO 95 J= IPs6 
IF(GOUDCJ) 295595292 
92 SUM = SUM*+P(LeTaJ) *COSFCWNCL) *#0C Is J) *COSFCTHETACK)@PSICI»J))) 
xX *PCLoJeT*SINFCWN(LI*DCTe J) *COSFCTHETACK )°PSICI9Jd))) 
95 CONTINUE 
90 CONTINUE 
60 ECLsK)=CAVEP(L)4+2e0#SUM)/TERMS 
IFCSENSE SWITCH 1)52963 
62 PRINTOsFQCL) oAVEPCLISCIDCI slEtoit ds CECLsK) sK31072) 
6 FORMATC AX 1PE11.4°1PEL2 4s 1 1A6/(1001X91PE11.4))) 
63 G=0.1/DTR 
DELTA THETAES5DEG UR H=1/2*DEL THETA 
DEC1)EGRCECL®2)°F(L»72)) 


DEC72I=G¥CE(CLo1 “ECL *71)) 


61 


66 


67 


69 


68 


64 


65 


DO 61 ITHETA22s71 
DECITHETAI=G*(ECLsTTHETA+1 ECL» ITHETA@1)) 
NEC73I=DEC1) 
K=4 

00 65 I=1972 


IF CDEC I) *DECI+1) 166267265 


FIV 1S 5 DEG IN RADYANS 
THZROCKI=FIVECZIFCDECID/C DECI )=DECI+1)))) 
K=K+]1 

GO TO 65 

IFCDECI) 16895956" 
THZROCK)=FIV*FLOATFCI“1) 
K=K+1 

IF CDEC I 41) 96526465 
THZROCK DEFIV*eFLOATFCI) 
K=K+]1 

CONTINUE 

NZEROS=K=1 

D0 75 K=ieNZEROS 

TERMS SAME AS ABOVE 
SUM2000 

DO 76 1=195 


IFCGOODCI) 176975977 


77 IP=Iel 
D0 79 J=IP»6 
IFCGQUD( I) 279979978 
78 SUM = SUMtPCLoT oJ) *#COSFCWNCL)*O0CIsJ)*COSFCTHZROCK)@PSICI»J))) 
x “P(LeJeTI*SINFCWNCL)*DCI 9d) *COSFCTHZROCKI@PSTC(I9Jd))) 
79 CONTINUE 
76 CONTINUE 
75 TEMCKI=CAVEP (CL) +2 09*SUMI/STERMS 
EMAX(L)=TEM(1) 
THMAKCLI=SRTU*CTWOPT@THZRO(1)) 
NO 74 K=2eNZEROS 
TFCEMAXCLISTEMCKI)I73974074 
73 EMAXC(L)=TEMCK) 
THMAX(L)=RTD*®CTWOPT=THZROCK)) 
THMAX JS BEARTNG FROM MAGNITIC NORTH 
74 CONTINUE 
SUM=0.0 
00 89 J=1%5 
IFCGIUDCI) 80789984 
81 IP=[+1 
DO 83 J=IPs6 
IFCGOUDCJ) 283983982 
82 SUM=SUM+P(LoToJ)#PCLoTo J) +PC Lodo eP(Lodel) 


83 CONTINUE 
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80 CONTINUE 
SSQSAVEPCL)*AVEPCLI*#2.0¥SUM 
H(L)2SSQ*TERMS#EMAX(L) *EMAXCL) 
ATILDASAVEP(LY/TWUPI 
HTILDASSSQ*ATILDA*ATILDA*TERMS 
H1(L)=1.0°CH(L)/HTILDA) 
H2CLYSCEMAXCLI=ATILDAJ/CAVEPCLI@ATILDA) 
IFCSENSE SWITCH 129970 
99 PRINT7>» TERMS*EMAX¢CLISATILDASHTILDA 
7 FORMATCIXsF5el1s1P3E1164 ) 
70 CONTINUE 
PRINT 35CID(I)»T=le11) »(GOUDCI)»T=196)0mM 
3 FORMATCIHL,39HLEAST SQUARE SINGLE WAVE TRAIN FIT OF ,11A6/ 
X1X96F2,092Xe2HM=I3// 
X 3X5 DHFREQUENCY 2 5X*S5HPERIOD» 3X91 1HNAVE LENGTH? 1Xs 1 LHATTENUATION? 4X» 
XSHAVE PsBXotHArlOXeSHBRNG »8Xa1HH912X »2HH1 oe L0X»2HH2) 
PRINT Gs CFQ(L)sPERTONCL) sWAVLGH(L?) »DEPATN(L) »AVEPCL) »EMAX(L)>» 
XTHMAXCL) HCL) HL CL) sH2CL) »LS25MP) 
Q@ FORMATCIX20PF 110679 3XS0PF9eSe1XSOPFItodsiXsiPEli ods 
X 1X9 1PE110491X91P E14 40 3X2 OP F729 3X0 1PE11 49 4Xs0PF6,356X20PF6.3 ) 
0G 110 L=2sMP 
AVEPCL) =E2TEN*LOGP CABSFCAVEPCL))) 
EMAX(L)SE2TEN*LOGF CABSFCEMAX(L))) 


110 WAVLGHCL)=EMAX(L)SE2TEN*LOGFCABSFCDEPATNCL))) 


AVEP(1)=AVEP(2) 
EMAXC1)=EMAX(€2) 
WAVLGHC1)=WAVLGHC 2) 
CALL SYMBOL(2.05%.75%0.1214HWAVE TRAIN.FN= 20.014) 
CALL NUMBER(3.5500752001%FQC(MP) 900092) 
CALL GPHPVW(MP» IN» AVEPs EMAX»WAVLGH) 
WAVLGH IS SURFACE WAVE STAFF SPECTRUM EST FROM BOTTOM PRESSUXE. 
AVE=P EMAy=V WAVLGH=w ON THE PLOTS, 
ZM=M 
NP=DELTAT*ZM*0.541.0 
DELX=0,.0 
90 115 I=156 
CALL NUMBER(DELX20.120.01»GO0D(1)»0.0s—1) 
115 DELX=DELX+0.1 
CALL PFB3DCNP* IDs EMAXsFQ*THMAX ) 
lei) fa) SS) 
100 REWIND 3 
CALL EXIT 
PAUSE 70707 
GO 79 15 
END(O%15190s1) 


END 
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FORTRAN TO ALGOL TRANSLATOR 
PHASE 1 FORTRAN STATEMENTS 


SUBROUTINE GPHPVW(L»ID»ZLOGP»ZLOGV.ZLOGW) 
SPECIAL FORM OF GPHPVW FOR 1335072 AUGUST 1967 


DIMENSION ZLOGPC2)57L0GVC2)»ZLOGW(2),10¢€2) 

CALL PLOTCO0.020.9%3) 

CALL PLOT (0.0219e592) 

CALL PLOT (8.091005s2) 

CALL PLOT (8.090.052) 

CALL PLOT (0.020,022) 

CALL PLOT (€2e0s1-59%3) 

CALL SYMBOL €0.0°71.097.1910(1)50.0%66) 

CALL AXIS (0.950.0,20HNORMALIZED FREQUENCY 5=20550002050+020.2) 
ZMAX2ZLOGP(1) 

DO 10 T=2°L 

COMPAR=ZLOGPCI) 

ZMAX=MAX1FCZM4X»CUMPAR) 

MAX=7MAX+3.0 

RE=SMAXK=8 

CALL AXISC0O+0%000247HLOG POWER DENSITY172820290¢0sBE5 100) 
DX=5,0/FLOATFCL=1) 

Y=ZLOGP(1)=BE 

X=0.0 


CALL PLOTCX»Y»3) 
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20 


21 


23 


DO 20 I[=2°L 

X=X+DX 

Y=ZLOGPC(I)=°B6E 

CALL PLOT (XsY»s2) 

CALL SYMBOL(XsY¥sOelelHPs0.001) 
X=5.0 

Y=ZLOGVC(LI=BE 

CALL SYMBOLCXsYeOolstHAs0.0s1) 
CALL PLOTCXsY»3) 

M=L=1 

DO 21 T=1»™ 

X=SKX=NX 

Itsiel 

YEZLOGVCII)@BE 

GALE PLOVGhovoe» 

WMAX=MAX 

X=0.0 

Y=ZLOGWC(1)°bE 

CALL PLOTCXsYe3) 

DO 22 T=2°L 

X=X+9X 

TF CZLUGWCI)@WMAX)24523923 
Y=8.0 


GO TO 22 
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24 
22 


30 


40 


50 


YeZLOGW(I)=BE 
CALL PLOT(CX»Y»2) 


CALL SYMBOL (Xs Y¥s%e151HS20,001) 


IF (SENSE LIGHT 


124030 


CALL PLOT (©2.099e0»923) 


SENSE LIGHT } 

GO TO 50 

CALL PLOT (6.02 
RETURN 

END( 991919029) 


END 


12.0973) 
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FORTRAN TO ALGOL TRANSLATOR 
PHASE 1 FORTRAN STATEMENTS 
SUBROUTINE PF33DCNrID»PyRoB) 
P=FUNCTION OF R = POWER DENSITY 
B= FUNCTION OF R = COMPASS BEARING 
R= FREQUENCY IN HZ 
DIMENSION J[DC2)2PC62)9RC2)5BC2) 966360) 95350) 
T=10.0 
D=0.70710676 
DTH=9 2017453293 
DTH ITS ONE DEGREEE IN RADIANS IE RADIANS PER 1 DEGREE 
A=0.9 
DO 10 1[=1»%360 
A=A+0TH 
CCI) =COSFCA) 
SCI)J=SSINFCA) 
CALL PLOTCO0.950.023) 
CALL PLOTC0,0%10.5»2) 
CALL PLOT(B.0910.592) 
CALL PLOT(8.050.922) 
CALL PLOTC0.0°0.092) 
CALL SYMBOL(1.09%.5570.121D(1)90.0266) 
CALL PLOTC4,053.52=°3) 


CALL PLOT(3.0°0.092) 
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30 


20 


CALL SYMBOL(3.25%0,020.121HN +0,021) 

CALL AX1ISC0.0°1.0°0H £095,0290.0%71.0%1.0) 
CALL PLOTC0.091.093) 

CALL PLOTC0.0»0.092) 

CALL PLOTC#3.020.092) 


CALL SYMBOL(23.2690,.0%0e1%1HS »0¢001) 


CALL SYMBOL(X*0.12Y"0els0e1s1HE » 00001) 
CALL PLOT(X»Ys3) 

CALL PLOTC#Xs=Y¥s>) 

CALL SYMBOL (O.1°XsM.1°Ys0e1ls1HW » 06081) 
CALL PLOT(0.020.02°3) 

DO 20 J=1%6 

A=0.5*FLOATFCI) 

CALL PLOTCA20.023) 

DO 30 J=i»360 

Y=A*S(CJ)*D 

XSA*CCJ) + Y 

CALL PLOTCXsYe2) 

CONTINUE 

CALL PLOTCO.020.99°3) 

D=10.0#D 


10*D0 NEEDED TO SCALE 0.070.3 TO 073 INCHES ON PLOTS 
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49 


42 


41 


50 


DO 40 J=2eN 
RAD=DTH*B(1) 
Y=RCT)*O*SINFCeORAD) 
XERCT)*T*COSFCRAD) + Y 
TUD=3 
CALL PLOTCXsYsIUD) 
AY=P(I1) + Y + 1.96 
2 IS NOT ADLED S9 ToP SYMBOL WILL BE THE REFERENCE IN 
CALL PLOTCXsAYs2) 
CALL SYMBOL (XsAYs0.082250,0%°2) 
CALL PLOTCXs AY» 3) 
CARE PEOn GX Yi92) 
IFCSENSE LIGHT 1941542 
CALL PLOTC£4.0>7.0,23) 
SENSE LIGHT 1 
GO 79 50 
CALL PLOTC4,05°14.0,°3) 
RETURN 


END( 0915912091) 
END 
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APPENDIX C 
A FORTRAN PROGRAM FOR ITERATIVE WAVE TRAIN ANALYSIS 


The FORTRAN listing of a Burroughs B-5500 program for the itera- 
tive least square multiple wave train analysis of a spectral matrix 
is presented. The mathematics is an iterative utilization of the 
single wave train analysis described in the body of this report. 
Following the single wave train analysis of the measured spectral 
matrix the resulting values of single wave train power, A, and wave 
bearing, 9, are used to find the spectral matrix that would occur for 
such a wave. Details for this are given in Section 5 of this report. 
From the above, a residual measured spectral matrix is formed by sub- 
tracting a fractional portion of the single wave spectral matrix from 
the previously used spectral matrix. The residual spectral matrix is 
then single wave train analyzed. The above procedure is continued 
iteratively until a specified number of iterations have been com- 
pleted or the residual spectral matrix total power gets smaller than 
a specified value. Table Cl shows numerical results for several fre- 
quencies. Figure Cl is a plot of the bearing, 6, and the ratio of A 
to the total power available in the original measure spectral matrix 
for the frequency band 0.00833 to 0.24187 Hz. The iteration parameters 
were set for a maximum of five iterations, a residual power ratio of 
Q.1, and a fractional portion value of 0.1. Both results are for data 
collected for task SWOC at Stage II between 1220 and 1249 hours on 
4 June 1965. The wind speed was 12 knots from a bearing of 280 degrees. 
In Table Cl A, H, and bearing are as previously defined. AVE P is the 
average C;,(f) for the residual spectral matrix and P is the AVE P for 
the first iteration; i.e., the original measured average power. The Hl 
and H2 values are measures of the isotropicity of the energy represented 
by the residual spectral matrix. Hl compares the least square value of H 
with the value H’ that would have been obtained if the wave energy were 


isotropic: 
ob eo ak cs (Geyfsll)) 


H2 compares the power A of a single wave fit to the total power, AVE P, 
and the value A’ that would have been obtained for isotropic energy: 


H2 = (A-A')/ (AVE P-A'). 


Both Hl and H2 are between 0 and 1, the lower limit is for the isotropic 
case and the upper for a plane wave from a single direction. 
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